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ABSTRACT 

Screening  methods  arc  beneficial  for  studies  involving  simulations  that  have  a  large  number  of  variables 
where  a  relatively  small  (but  unknown)  subset  is  important.  In  this  paper,  we  show  how  a  newly  proposed 
Lasso-optimal  screening  design  and  analysis  method  can  be  useful  for  efficiently  conducting  simulation 
screening  experiments.  Our  approach  uses  new  criteria  for  generating  supersaturated  designs,  and  a  new 
algorithm  for  selecting  the  optimal  tuning  parameters  for  Lasso  model  selection.  We  generate  a  24x69 
Lasso  optimal  supersaturated  design,  illustrate  its  potential  with  a  numerical  evaluation,  and  apply  it  to  an 
agent-based  simulation  of  maritime  escort  operations  in  the  Strait  of  Gibraltar.  This  application  is  paid  of  a 
larger  project  that  seeks  to  leverage  simulation  models  during  the  ship  design  process,  and  so  construct  ships 
that  are  both  cost  effective  and  operationally  effective.  The  supersaturated  screening  design  has  already 
proved  beneficial  for  model  verification  and  validation. 

1  INTRODUCTION 

Complex  systems  such  as  large-scale  computer  simulation  models  typically  involve  a  large  number  of 
factors.  When  investigating  such  systems,  screening  experiments  are  often  used  to  sift  through  these  factors 
and  identify  a  subgroup  that  most  significantly  influences  the  response  of  interest.  A  secondary  experiment 
can  efficiently  allocate  the  remaining  computational  resources  to  those  important  factors,  and  this  two-phase 
process  can  therefore  greatly  expand  the  ability  of  analysts  and  decision  makers  to  gain  insights  about  a 
complicated  system  with  a  reasonable  computing  budget. 

The  focus  of  the  paper  is  on  using  supersaturated  designs  for  factor  screening.  A  supersaturated  design 
is  a  fractional  factorial  design  with  more  factors  than  experimental  runs  (Booth  and  Cox  1962,  Lin  1993, 
Wu  1993).  Supersaturated  designs  share  the  same  basic  principle  (i.e.,  the  sparsity  of  effects  principle) 
as  screening  experiments,  and  arc  specifically  attractive  for  high  dimensional  scenarios  since  they  require 
far  less  runs  than  the  number  of  factors.  On  the  other  hand,  supersaturated  designs  do  not  have  sufficient 
degrees  of  freedom  to  simultaneously  estimate  all  main  effects  for  all  factors.  Therefore,  different  analysis 
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methods  arc  needed  than  those  used  for  non-saturated  designs.  Recently,  L\ -penalty-based  methods  such  as 
Lasso  have  been  developed  and  shown  great  potential  for  analyzing  high  dimensional  data  (Osborne  et  al. 
2000,  Efron  et  al.  2004).  Our  latest  work  (Xing  et  al.  2013a, b)  discusses  a  new  approach  to  utilize  Lasso  for 
factor  screening  and  proposes  new  searching  criteria  to  construct  Lasso-optimal  supersaturated  design.  In 
this  paper,  we  review  the  method  and  generate  an  optimal  supersaturated  design  to  explore  an  agent-based 
simulation  model  of  marine  escort  operations.  The  agent-based  simulation  model  involves  a  naval  escort 
mission  scenario  in  an  anti-surface  warfare  environment.  It  is  based  on  a  2002  Morocco  incident,  where 
terrorists  planned  to  conduct  suicide  attacks  on  vessels  passing  through  the  Strait  of  Gibraltar.  The  screening 
experiments  take  24  runs  to  study  65  factors,  and  our  proposed  Lasso  method  picks  up  to  11  important 
factors. 

The  paper  is  organized  as  follows.  Section  2  focuses  on  the  methodology  used  in  this  paper.  We  review 
the  Lasso  method  and  its  application  in  factor  screening,  discuss  the  probability  of  correct  model  selection 
using  Lasso,  propose  a  new  criterion  to  generate  Lasso-optimal  supersaturated  designs,  and  discuss  the 
determination  of  the  Lasso  tuning  parameter.  In  Section  3,  we  describe  a  simulation  model  of  maritime 
escort  operations,  and  study  the  model  with  a  24  by  65  Lasso-optimal  supersaturated  design.  We  conclude 
with  a  brief  discussion. 


2  METHODOLOGY 

We  begin  by  briefly  describing  the  Lasso-optimal  screening  design  and  analysis  method  used  in  the  two 
paper  to  analyze  the  simulation  model.  Lor  more  details,  please  refer  to  Xing  et  al.  (2013b)  and  Xing  et  al. 
(2013a). 


2.1  LASSO  INTRODUCTION 

The  least  absolute  shrinkage  and  selection  operator,  or  Lasso  in  short,  refers  to  a  group  of  methods  that  use 
an  L\  penalty  to  shrink  parameter  estimates  and  perform  automatic  variable  selection  (Tibshirani  1996). 
Lasso  does  both  continuous  shrinkage  and  automatic  variable  selection  simultaneously,  and  has  shown 
great  potential  for  analyzing  high-dimensional  data. 

Consider  a  linear  regression  model 

Y  =  X/3  +  £, 

where  Y  =  (Y\,Y2,  ■ .  . ,  Yn)'  is  the  n  dimensional  vector  of  observed  responses,  X  is  the  n  x  p  design  matrix, 
jS  =  ()3i ,  /3t  ,  -  - . ,  ^p/isap-dimensional  vector  of  the  unknown  regression  coefficients,  ande  =  (£1,62,  •  •  ■ ,  £n)' 
is  the  vector  of  random  errors  from  a  multivariate  normal  distribution  N(0,<J2In).  Denote  the  columns 
of  X  as  (XhX2,...,Xp);  these  represent  the  p  independent  variables.  Lor  i  =  1 . 2. .../;,  X,  is  said  to  be  a 
true  (important)  variable  if  /J,  /  0;  otherwise,  X\  is  said  to  be  an  untrue  (unimportant)  variable.  The  Lasso 
estimate  of  /i,  denoted  by  /m,  is  the  solution  of  the  following  L\  -penalized  least  squares  problem: 


5A  A 

p  =  -  arg  mm 

2  BeRp 


£  L-LS/Xy  +  A£|  Bt 

1=1  \  7=1  /  i=l 


(1) 


where  |B,|  is  the  Li-norm  penalty  function  of  B,  and  X  is  the  tuning  parameter  that  controls  the  amount 
of  penalty.  Note  that  we  can  assume,  without  loss  of  generality  for  factor  screening,  that  the  T’s  have  been 
centered  with  mean  zero  so  the  intercept  is  also  equal  to  zero  (Tibshirani  1996).  Lor  ease  in  discussion,  let 
A  =  {/ :  /j,  /  0}  and  C  =  {i  :  /3,  =  0}  be  the  collection  of  the  true  and  untrue  variables  indexes,  respectively. 

Given  a  specific  A  >  0,  there  are  many  fast  algorithms  to  solve  Lasso,  such  as  LARS  (Efron  et  al. 
2004),  Coordinate  Descent  (Lriedman  et  al.  2007),  and  others.  In  this  paper,  we  use  the  LARS  algorithm 
because  of  its  demonstrated  efficiency  for  large-scale  data.  Lor  a  fixed  A,  define  AA  =  {  j  :  (i2  /  0},  which 
is  referred  as  the  collection  of  selected  variables  (those  estimated  to  be  important  factors)  at  A.  Denote 
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A  Q  A  T  ^  1 

C  as  the  complement  of  A  .  When  A  increases,  more  p;-  s  will  become  zero  due  to  the  heavier  impact 
of  the  penalty  term  in  Equation  (1).  The  path  of  fVL  as  A  changes  is  referred  to  as  the  solution  path  of 
Lasso  (Tibshirani  1996),  or  simply  the  Lasso  path.  Different  methods  have  previously  been  proposed  for 
choosing  the  appropriate  tuning  parameter  X  which,  in  turn,  controls  how  many  variables  will  be  selected. 


2.2  PROBABILITY  OF  SIGN-CONSISTENT  SELECTION  VIA  LASSO  AND  LASSO-OPTIMAL 

SUPERSATURATED  DESIGNS 

For  any  A  >  0.  /J  'k  will  always  be  a  biased  estimator  of  the  true  parameter  vector  /3.  There  are  two  ways  to 
evaluate  the  Lasso  estimate  j8a.  One  is  to  see  whether  the  non-zero  items  in  match  the  true  important 
factors,  which  we  call  variable  selection  consistency;  and  another  is  to  measure  the  difference  between  the 
estimated  coefficients  and  /J.  Here  we  focus  on  the  first  one,  which  is  consistent  with  our  factor  screening 
goal. 

There  arc  two  types  of  variable  selection  consistency  for  Lasso  discussed  in  the  literature.  When 
A  =  A,  that  is,  all  the  true  variables  arc  selected.  Lasso  is  said  to  have  ordinary  consistency  in  variable 
selection  at  A.  Similarly,  Lasso  is  said  to  have  sign  consistency  in  variable  selection  at  A  if  the  p  x  1 
vectors  of  true  and  estimated  signs  of  the  coefficients  arc  equal.  Mathematically,  define  the  function  sgn(B) 
to  operate  element-by-element  on  vector  B  so  that  sgn(Bj)  =  takes  on  values  in  {—1,0, 1}.  Then 

sign  consistency  means  that  sgn( /D)  =  sgn ( /j  ) .  Notice  that  the  sign  consistency  of  Lasso  implies  ordinary 
consistency,  but  the  converse  is  not  generally  true. 

Let  q  denote  the  number  of  true  variables.  Let  Pa  be  the  subvector  of  p  that  corresponds  to  the  true 
variables,  let  Xa  be  the  submatrix  of  X  with  columns  in  A,  and  let  Xc  be  the  complement  of  Xa-  From 
these,  define  P  =  XA(XAXA)~xXA,  R  =  X^XA(XAXA)~xsgr\(PA),  and  D  =  ^(XAXAylsgn{PA)- 

The  following  theorem  by  Xing  et  al.  (2013b)  gives  the  probability  that  Lasso  achieves  sign  consistency 
for  a  given  sample  of  size  n,  design  matrix  X,  and  penalty  parameter  A. 

Theorem  1.  Let  W  be  a  (p-q)-dimensional  multivariate  N(RAX'C(I  —  P)Xc<J2/X2),  and  let  V  be  a  q- 
dimensional  multivariate  N(D ,  (XAXA)~  1  c2).  The  probability  that  Lasso  achieves  sign  consistent  selection, 
which  is  denoted  by  P(SCS),  is: 

P(SCS)  =  P(<fi)P(<f2),  (2) 


where  S\  and  Si  are  the  events 


Si  ={-l  <W<  1}  and 


Si  =  f]  {Pi  i  [min(0,  V;),  max(0,  V,)]} . 

ieA 


With  fixed  A,  /3,  A  and  a2,  the  probabilities  of  S\  and  Si  depend  on  the  following  four  terms:  R  and 
D ,  which  are  the  mean  vectors  of  IT  and  V,  along  with  4X^(7  —  P)Xc<J2  /  X2  and  ^(XAXA)~  'a2,  which  are 
the  variances  of  W  and  V,  respectively.  If  we  can  optimize  P(SCS),  we  will  optimize  the  performance  of 
Lasso  for  factor  screening.  There  are  two  challenges,  however.  First,  the  direct  evaluation  and  optimization 
of  P(S\ )  and  P (Si)  with  respect  to  P ,  A,  and  X  require  high-dimensional  integration,  which  becomes 
computational  prohibitive  as  the  dimensions  increase.  Second,  these  probabilities  depend  on  unknown 
parameters. 

To  address  the  first  challenge,  Xing  et  al.  (2013b)  proposed  both  a  lower  bound  method  and  Monte  Carlo 
simulation  method  to  approximate  P(S] )  and  P {Si)  efficiently  in  high  dimensional  cases.  Notationally,  let 
Dj  {i  =  1, . . . ,q)  denote  the  7th  element  of  D ,  and  let  Rj  —  q)  denote  the  jth  element  of  R.  Let 

( ,  (J2y2 , . . . ,  Cjp  J  be  the  vector  of  diagonal  elements  of  (4X'c(l  —  P)Xc<72),  and  let  (o^ ,  ap2 , . . . ,  op  / 
be  the  vector  of  the  diagonal  elements  of  ({X'aXa)~{  o2).  Finally,  let  /(•)  be  an  indicator  such  that  I (x)  =  1 
if  x  is  true,  and  I(x)  =  0,  otherwise.  The  lower  bound  method  is  based  on  the  following  proposition  from 
Xing  et  al.  (2013b): 
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Proposition  1.  A  lower  bound  for  P(S\ )  is 


P(^i)w  =  max 


0,1 "  <  >  ^C/)) 


1 


(3) 


where 


1 


Ql{i)=U+Rj)eXp' 


-Z2(l+Rj 

2c7w. 


+  (T ■=*j)‘qr 


-A2(l-Rj 

2<7Wt 


and 


QiU)  = 


-exp 


-Z2(\+Rj 


V2n(l  +Rj) 

Also,  a  lower  bound  for  is 


2<7W; 


1 

+  I“' 


v^f  (1  -Rj)2  +  if 


exp  ■ 


-A2(l  —Rj 

2(%, 


P(A)loW  =  max  [0,1-  oVi  ^  ^ 


i—l 


l 


I  ft  I  V  2ft -AD, 


exp 


~(2ft -AP, 
8  a2. 


(4) 


77?e«,  a  lower  bound  for  P(SCS)  is  given  by  P(SCS);ow  =  P(<^i)/0wP(<^2)/ow 

The  lower  bound  approximation  above  can  be  calculated  very  efficiently  compared  to  P(SCS).  The 
simulation  results  show  that  when  the  total  number  of  variables  p  <  25,  the  lower  bound  offers  a  satisfactory 
approximation  for  the  probability  of  correct  selection.  For  higher  order  cases,  however,  the  lower  bound 
deviates  substantially  from  the  true  value  and  the  approximation  of  P(Z'i )  and  P(<^)  may  approach  zero. 
In  these  cases,  we  use  a  straightforward  Monte  Carlo  sampling  method  instead  of  the  lower  bounds  in 
(3)  and  (4).  In  summary,  the  algorithm  sequentially  and  independently  generates  sample  vectors  from  W 
and  V,  both  of  which  are  multivariate  normal  random  variables,  to  obtain  direct  estimates  P(<P])  Monte  and 
P(^2 ) Monte,  respectively.  The  algorithm  stops  when  the  estimation  converges.  For  more  details,  please  refer 
to  Xing  et  al.  (2013b). 

The  second  challenge  of  directly  optimizing  P(SCS)  is  that  the  probability  depends  on  unknown 
parameters,  including  the  number  of  true  factors,  the  magnitudes  and  signs  of  [X\ ,  and  the  variance  a2. 
These  parameter  values  arc  not  available  in  the  screening  stage — if  they  were,  there  would  be  no  need  for 
a  screening  experiment.  One  approach  to  solve  this  problem  is  to  focus  only  on  the  influence  of  the  design 
structure  to  the  probability  of  correct  selection.  Specifically,  Zhao  and  Yu  (2006)  prove  that  if  \R\  <  1, 
then  P(SCS)  =  1  asymptotically.  This  is  called  the  irrepresentable  condition.  From  the  P(SCS)  expression 
in  Theorem  1,  it  is  easy  to  see  that  R  dominates  P(Z) ).  When  n  goes  to  infinity,  both  P(S\ )  and  P(Sf) 
will  converge  to  1  (independent  of  \R\)  if  A'j  <  I  ,  which  implies  the  irrepresentable  condition  discussed 
in  Zhao  and  Yu  (2006).  This  motivates  us  to  propose  a  standard  for  optimal  design  that  is  based  on  the 
residual  sum  of  squares  RSS  as  a  function  of  qmax,  the  maximum  number  of  potentially  important  factors 
defined  by  the  user.  For  ease  of  notation  we  will  drop  the  subscript  and  simply  use  q.  Let  A  be  a  submatrix 
of  selected  variables  (important  factors)  at  A,  and  define 


RSS  (A)  =  trace{X'cXA{X'AXA)-2X'AXc). 

Since  we  do  not  have  prior  knowledge  of  which  q  factors  are  important,  all  possible  subsets  of  q 
factors  need  to  be  considered  when  constructing  the  optimal  design.  Therefore,  we  evaluate  the  suitability 
of  design  X  for  selecting  q  factors  as  follows:  we  compute 

RSS,(X)  =  fP)  £  RSS  (A), 

W  a£/ 
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where  sZ  is  the  collection  of  all  possible  subsets  of  q  important  factors  among  the  p  factors.  For  fixed  p 
and  n,  we  refer  to  a  design  that  minimizes  RSS^(X)  over  all  possible  designs  with  specified  q  as  an  optimal 
RSS?  design. 

The  design  is  generated  by  a  partial  gradient  column  exchange  heuristic  algorithm  discussed  in  Xing 
et  al.  (2013a).  The  intuition  of  the  algorithm  is  to  find  the  partial  gradient  of  the  objective  function  with 
respect  to  each  column  in  X,  and  use  the  gradient  information  to  guide  the  column  exchange  by  moving 
the  least  favorable  column  in  the  desired  direction.  When  the  progress  stagnates,  fresh  columns  will  be 
introduced.  Due  to  the  limitations  of  the  space,  we  omit  the  details  of  the  algorithm  here  and  direct  the 
interested  readers  to  Xing,  Wan,  and  Zhu  (2013a).  We  have  run  extensive  numerical  evaluation  to  compare 
the  RSS-optimal  supersaturated  designs  with  other  supersaturated  designs  proposed  in  the  literature.  The 
RSS-optimal  designs  do  have  higher  probability  of  correct  model  selection  compared  to  other  designs  when 
Lasso  is  used  as  the  analysis  method. 

2.3  MODEL  SELECTION  USING  LASSO:  DETERMINE  THE  TUNING  PARAMETER  A 

With  the  optimal  design,  we  collect  observations  from  simulation  model  and  apply  the  LARS  algorithm 
proposed  by  (Efron  et  al.  2004)  to  analyze  the  data.  One  issue  that  has  not  yet  been  discussed  in  this  paper 
is  the  selection  of  A,  the  tuning  parameter  of  Lasso.  This  matters  because  the  magnitude  of  A  determines 
how  many  coefficients  will  be  non-zero,  therefore  the  amount  of  screening  that  will  occur.  We  recently 
proposed  an  algorithm  based  on  self-voting  using  ordinary  least  squares  estimation,  called  SVOE,  to  select 
the  optimal  A.  We  briefly  summarize  it  here.  For  more  details,  please  refer  to  Xing  et  al.  (2013b). 

To  gain  insight  into  the  relationship  between  A  and  j3  ,  consider  the  following  generic  example  with 
16  variables.  The  underlying  true  model  is  Y  =X/3  +£,  where  /3  =  (2, 0,2, 0,0, 0,0, 0,2, 0,0, 0,0, 0,0,0)'. 
Note  that  the  true  (important)  variables  in  this  model  are  factors  1,  3  and  9.  Our  X  is  the  12  x  16 
supersaturated  design  given  in  Table  4  of  Li  and  Wu  (1997),  and  we  let  £  ~  N(0,  cr2/^).  We  generate  one 
response  independently  for  each  design  point  using  Monte  Carlo  simulation  and  apply  the  LARS  algorithm 
repetitively  with  A  changing  from  70  to  0.  This  gives  us  the  Lasso  solution  path  of  Figure  1.  Here  the 
x-axis  is  the  A;  and  the  y-axis  is  the  estimated  /3  value.  When  A  >  70,  all  estimated  coefficients  are  0.  As 
A  decreases,  factors  are  added  to  the  model  one  by  one. 


Figure  1:  Lasso  path  for  a  generic  12  x  16  supersaturated  design. 

The  seven  lines  in  Figure  1  represent  the  sequential  changes  of  /3,  for  the  first  seven  factors  selected 
into  the  model;  we  stop  at  seven  factors  to  avoid  crowding  the  figure,  but  twelve  factors  are  selected  at 
A  =  0.  The  vertical  lines  split  the  figure  into  different  regions,  and  the  numbers  in  squares  specify  how 
many  factors  arc  selected  within  the  regions.  For  example,  in  Region  3,  where  A  is  between  10.7  and  16.5, 
three  factors  are  included  in  the  model  (factors  1,  3,  and  9).  Thus,  for  each  A  value,  there  is  a  corresponding 
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if  we  treat  this  as  if  it  is  the  true  /3  and  plug  it  into  Equation  (2)  in  Theorem  1,  we  can  solve  for  a 
corresponding  Xopl .  We  can  repeat  this  process  until  the  X  and  its  counterpart  Xopt  match,  and  treat  this 
converged  X  as  the  true  Xopt . 

We  call  our  heuristic  SVOE,  for  self-voting  based  on  OLS  estimation.  It  is  based  on  the  self-voting 
principle  discussed  in  Gu  (1992).  The  numerical  evaluation  in  Xing  et  al.  (2013b)  shows  that  this  heuristic 
works  very  well  except  in  some  rare  cases.  When  there  arc  multiple  self-voting  Xopt,  the  user  can  determine 
which  model  they  would  like  to  pick.  Compared  with  other  parameter  tuning  methods,  for  example,  CV 
(Devijver  and  Kittler  1982),  AIC  (Akaike  1974),  BIC  (Schwarz  1978),  that  focus  on  finding  the  balance 
between  goodness  of  fit  and  overfitting  by  adding  penalties  for  selecting  too  many  factors  into  the  model, 
our  approach  focuses  on  directly  optimizing  the  probability  of  the  selecting  the  true  model.  Numerical 
evaluation  shows  that  our  SVOE  method  selects  the  best  X  the  majority  of  the  time  (Xing  et  al.  2013b), 
and  we  apply  it  in  the  next  section. 

3  MARITIME  ESCORT  OPERATIONS  SIMULATION 

To  illustrate  the  potential  benefits  of  this  screening  procedure,  we  apply  it  to  study  an  agent-based  model 
of  maritime  escort  operations.  A  brief  summary  of  the  context  follows.  For  further  details,  see  Kaymal 
(2013):  his  thesis  is  paid  of  a  multinational,  multi-year  effort  to  leverage  modeling  and  simulation  in  order 
to  better  understand  trade-offs  between  cost  and  operational  effectiveness  for  surface  ships,  and  so  improve 
the  ship  design  process. 

3.1  MOTIVATION 

In  the  past  two  decades,  there  has  been  a  significant  shift  in  navy  missions  towards  operations  other  than 
war.  Counter-piracy,  search  and  rescue,  maritime  interdiction,  maritime  patrol,  naval  escort  operations,  and 
humanitarian  efforts  are  the  main  focus  of  most  fleets  today,  but  the  vessels  that  are  currently  being  used  in 
such  operations  were  mainly  built  for  other  purposes.  For  instance,  on  17  August  2009,  the  North  Atlantic 
Council  approved  ’’Operation  Ocean  Shield”  to  fight  piracy  in  the  Gulf  of  Aden.  Among  six  surface  ships 
that  were  assigned  in  January-June  2012  rotation  of  this  NATO  mission,  one  was  a  destroyer  and  three 
were  frigates.  These  arc  sophisticated  warships  which  arc  capable  of  anti-surface  warfare,  anti-air  warfare, 
and  anti-submarine  warfare.  Although  these  sophisticated  multi-mission  capable  fleets  are  able  to  achieve 
good  results  in  expeditionary  warfare  against  a  strong  enemy  (Murphy  2007),  their  full  capabilities  will 
probably  be  used  in  less  than  1%  of  their  total  life  time.  Frigates  and  destroyers  arc  also  expensive  to  build 
and  operate. 

Smaller  combatants  are  much  cheaper  and  better  suited  for  the  more  common  modern  naval  operations 
due  to  their  flexibility.  Therefore,  many  nations  arc  reshaping  their  fleets  to  meet  emerging  operational 
demands,  by  stalling  to  build  multi-mission  capable  combatants  such  as  Offshore  Patrol  Vessels  (OPVs) 
that  arc  smaller  and  less  sophisticated  than  frigates  and  destroyers.  These  nations  arc  also  developing  new 
tactics  and  counter  measures  to  better  deal  with  the  new  threats.  The  main  purpose  of  OPVs  is  maintaining 
maritime  security  inshore  and  offshore,  and  they  can  be  deployed  globally  (Kimber  and  Booth  2010). 
Being  cheaper  yet  flexible,  OPVs  arc  also  great  options  for  those  countries  that  either  do  not  really  need 
or  cannot  afford  sophisticated  naval  combatants. 

Unfortunately,  the  ship  design  and  the  acquisition  process  has  not  kept  pace  with  the  rapid  technological 
advances  of  the  past  few  decades  (Ryan  and  Jons  2004).  Historical  cost  estimation  models,  such  as  those 
based  on  cost  per  ton,  often  fail  to  adequately  capture  new  technology-driven  requirements.  In  the  past, 
operational  requirements  were  often  not  considered  until  late  in  the  ship  design  process,  and  then  might 
be  based  on  limited  input  by  small  panels  of  subject  matter  experts.  Cost  effectiveness  and  operational 
effectiveness  arc  both  important,  and  it  is  extremely  hard  to  achieve  both  when  using  a  traditional  ship 
design  process.  Advances  in  modeling  and  simulation,  engineering  design  software,  and  massive  amounts 
of  data  now  allow  for  a  much  broader  and  richer  range  of  ship  specifications  and  capabilities  to  be  explored. 
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Moreover,  utilizing  simulation  and  analytical  models  to  build  decision  making  tools  will  ensure  collaboration 
between  warfighters  and  engineers  in  the  early  stages  of  the  process.  Exploiting  technology  is  paramount 
for  accomplishing  the  navy  objectives  and  increasing  the  effectiveness  for  both  cost  and  operations  (Mizine 
et  al.  2012). 

3.2  SIMULATION  SCENARIO  DESCTIPTION 

Even  though  the  term  ’’maritime  terrorism”  has  been  in  the  literature  for  more  than  a  few  decades,  it  was 
not  really  spelled  out  clearly  until  a  series  of  incidents  started  in  2000.  In  January  2000,  terrorists  tried  to 
attack  USS  The  Sullivans  (DDG-68)  in  Yemen.  They  failed  because  the  terrorist  boat  sank  right  before 
the  attack.  Following  this  event,  the  terrorists  attacked  USS  Cole  (DDG-67)  with  a  suicide  boat,  killing 
17  of  the  crew  members,  in  October  2000.  Almost  two  years  after  these  incidents,  in  October  2002,  a 
boat  with  explosives  hit  the  French  oil  tanker  Limburg  which  was  close  to  Yemen  coastal  waters  (Luft  and 
Korin  2004).  These  are  just  a  few  examples  of  how,  with  just  a  small  boat,  terrorists  can  cause  damage 
to  multi-million/billion  dollar  ships,  and — more  importantly — kill  innocent  people.  Many  other  terrorist 
activities  which  were  discovered  in  their  planning  phases  and  prevented.  These  incidents  show  that  maritime 
terrorism  is  a  serious  problem,  and  that  surface  combatants  must  be  ready  to  fight  terrorists  at  sea. 

In  June  2002  a  group  of  terrorists,  who  were  planning  an  attack  on  two  merchant  vessels  in  the  Strait 
of  Gibraltar,  were  caught  by  Moroccan  officials  (Maggio  2008).  The  fact  that  the  officials  arrested  those 
terrorists  does  not  mean  that  similar  plans  will  not  be  put  into  practice  by  other  groups.  If  a  terrorist  group 
fills  a  boat  with  explosives  and  approach  a  ship  in  Strait  of  Gibraltar-,  it  might  be  really  hard  to  identify 
that  boat  as  a  terrorist  boat  since  many  vessels  are  in  the  strait  at  any  point  in  time. 

Our  specific  scenario  involves  a  naval  escort  mission  scenario  in  an  anti-surface  warfare  environment. 
It  is  based  on  the  2002  Morocco  incident,  and  implemented  in  an  agent-based  simulation  modeling  platform 
called  MANA,  designed  by  the  Operational  Analysis  personnel  of  the  New  Zealand  Defence  Technology 
Agency  (DTA)  (Lauren  and  Stephen  2002).  The  key  attributes  of  MANA  which  make  it  a  useful  tool 
for  military  applications  are  the  situational  awareness  of  the  agents,  advanced  communication  capabilities 
within  squads  and  with  other  agents,  the  interaction  of  entities  with  friends  and  foes,  and  the  user  friendly 
design  of  the  program.  There  are  four  main  sets  of  parameters  that  form  agent  behaviors  in  MANA: 

•  Personality  weightings  determine  the  willingness  of  agents  to  perform  a  particular  action; 

•  Movement  constraints  modify  the  basic  personality  weightings  of  the  agents; 

•  Intrinsic  capabilities  determine  the  physical  characteristics  of  the  agents  such  as  sensors,  weapons, 
or  fuel  level;  and 

•  Movement  characteristic  adjustments  ensure  that  agent  actions  change  in  different  terrain  conditions 
and  different  situations. 

Just  recently,  DTA  released  MANA-V.  The  “V”  stands  for  both  vector  and  five.  In  this  version,  the 
programmers  replaced  the  cell-based  movement  of  previous  versions  with  vector  based  movement.  This 
allows  building  larger  battlefield  regions  with  panning  and  zooming  options,  and  allows  defining  users  to 
define  distances,  and  the  attributes  such  as  speed  and  range,  in  real-world  units  rather  than  pixel-based 
units  (McIntosh  2009). 

There  are  six  types  of  agents  in  the  scenario:  the  high  value  unit  (HVU)  being  escorted,  the  OPV,  the 
OPV’s  helicopter,  terrorist  boats,  known  vessels,  and  unknown  vessels.  A  screenshot  from  MANA  appeal's 
in  Figure  1;  note  that  the  sizes  of  the  agents  are  not  to  scale,  but  are  magnified  for  easier  visualization  by 
the  user.  The  allegiance  of  the  HVU,  OPV,  and  the  helicopter  is  “friend,”  the  allegiance  of  the  terrorist 
boats  is  “hostile,”  and  the  allegiance  of  the  known  and  unknown  vessels  is  “neutral.”  The  OPV’s  mission 
is  to  escort  the  HVU  through  the  Strait  of  Gibraltar-  and  protect  the  HVU  from  attacks  that  can  occur  in 
the  passage.  The  OPV  has  several  guns,  ranging  from  high  caliber  guns  to  machine  guns.  It  also  has  a 
helicopter  landing  platform.  The  main  role  of  the  helicopter  is  to  detect  and  classify  unknown  vessels.  Its 
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high  speed  and  maneuverability  give  the  friendly  forces  an  advantage  against  the  hostiles.  The  helicopter 
can  also  be  equipped  with  a  machine  gun  so  that,  if  necessary,  it  can  start  firing  before  the  hostile  vessels 
come  too  close  to  the  HVU. 


Figure  2:  MANA  screenshot,  with  six  types  of  agents  in  the  Strait  of  Gibraltar  (from  Kaymal  2013). 

The  terrorist  boats  arc  loaded  with  explosives,  and  their  purpose  is  to  get  close  enough  to  the  HVU  to 
perform  a  suicide  attack.  The  only  target  of  the  hostile  boats  is  the  HVU.  They  do  not  attack  the  OPV,  but 
in  some  cases  they  may  try  to  evade  it.  One  of  the  critical  properties  of  the  terrorist  boats  is  that  they  are 
initially  acquired  as  unknown  vessels  in  the  friendly  force  radar  systems  until  they  are  classified  as  enemy. 
Therefore,  either  the  helicopter  classifies  them  as  enemy,  or  they  get  closer  to  the  OPV  or  HVU,  and  the 
friendly  ships  classify  them  depending  on  the  range. 

The  known  ships  and  the  unknown  ships  are  both  neutral  and  they  do  not  pose  a  threat  to  friendly 
assets.  The  only  difference  between  the  two  is  that  while  the  OPV  can  instantly  classify  the  known  neutral 
ships  using  an  Automatic  Identification  System  (AIS)  device,  unknown  ships  cannot  be  classified  when 
they  are  initially  detected.  This  occurs  for  several  reasons:  some  ships  may  be  too  small  to  carry  an  AIS 
device,  others  may  have  a  device  that  is  not  working  or  turned  off.  This  represents  the  type  of  marine 
traffic  often  found  in  the  Strait  of  Gibraltar  and  other  similar  situations.  It  complicates  matters  for  the 
OPV,  since  the  OPV  and  its  helicopter  may  be  unable  to  fire  on  hostile  boats  if  neutral  ships  are  close  by. 
Moreover,  the  OPV  and  helicopter  may  need  to  spend  a  substantial  amount  of  time  and  effort  identifying 
unknown  ships,  even  though  the  vast  majority  are  neutral.  These  distractions  might  give  terrorists  a  better 
chance  to  approach  the  HVU  without  being  classified  as  threats. 

The  simulation  begins  with  the  HVU  and  OPV  approaching  the  strait.  It  halts  the  instant  that  either 
the  HVU  successfully  navigates  the  strait,  or  it  is  attacked  by  a  hostile  ship. 

3.3  SIMULATION  EXPERIMENTS 

We  conducted  multiple  sets  of  experiments  in  the  process  of  creating  and  exploring  this  simulation  model. 
Those  for  exploring  the  trade-offs  used  nearly  orthogonal  and  balanced  mixed  designs  (Vieira  Jr  et  al.  2011, 
Vieira  Jr.  et  al.  2013),  which  are  space-filling  designs  capable  of  handling  a  mix  of  continuous,  discrete, 
and  qualitative  factors.  The  initial  experiment  involved  100  replications  of  a  NOAB  design  for  38  factors, 
including  one  that  represented  the  time  step  used  in  the  scenario.  This  directly  affects  the  time  it  takes  to 
complete  the  simulation  runs,  but  the  time  step  must  be  small  to  avoid  model  artifacts. 

Once  the  choice  of  time  step  was  settled  and  some  other  minor  changes  were  made,  the  final  space-filling 
experiment  involved  200  replications  of  a  NOAB  for  35  factors.  Consolidating  these  results  into  a  single 
file  for  analysis  resulted  in  153,600  simulation  runs.  Kaymal  (2013)  uses  a  range  of  statistical  techniques, 
including  regression,  logistic  regression,  and  partition  trees,  to  model  the  relationship  between  the  factors 
and  key  responses.  These  models  also  are  used  to  suggest  factor  ranges  and  combinations  that  make  efficient 
use  of  resources  but  achieve  high  operational  effectiveness. 
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In  this  paper,  we  focus  on  just  one  of  the  many  potential  output  measures.  Our  response  Y  is  the 
proportion  of  time  that  the  HVU  successfully  navigates  the  strait.  We  use  a  24  x  69  Lasso-optimal  design 
capable  of  examining  up  to  69  factors,  and  identifying  up  to  24  important  ones.  The  design  is  available 
on  request.  We  augmented  KaymaTs  list  of  35  factors  with  an  additional  30  that  had  not  initially  been 
explored.  To  maintain  the  spirit  of  a  screening  experiment,  we  made  very  tiny  changes  in  factor  settings  for 
personality  weights  and  communication  characteristics  (e.g.,  ±2  on  a  scale  of  —100  to  100).  In  addition 
to  serving  as  an  example  for  the  Lasso  screening  approach,  we  viewed  this  as  an  opportunity  to  beta-test 
some  aspects  of  MANA  that  were  not  used  in  Kaymal’s  experiments.  (Over  the  years,  the  developers  have 
added  many  features  to  MANA  in  response  to  requests  by  NPS  students.  In  turn,  the  beta-testing  by  NPS 
students  of  MANA’s  features  has  been  useful  for  verification  and  validation  efforts.) 

To  test  the  performance  of  the  selected  design,  we  ran  a  quick  numerical  evaluation  of  the  selected 
design.  Consider  a  linear  model  with  65  factors  and  g  =  5, 10, 15,  and  20  respectively.  The  important  factors 
are  randomly  selected  and  their  effects  are  set  at  3.  The  other  factors  all  have  an  effect  of  0.  The  white 
noise  follows  a  standard  normal  distribution.  We  generate  one  set  of  responses  of  the  24  x  69  Lasso-optimal 
design  for  each  case,  and  the  results  arc  summarized  in  Table  1.  All  factors  that  arc  both  truly  important 
and  selected  by  Lasso  arc  shown  in  boldface.  Note  that  for  the  q  =  10  case,  there  arc  multiple  optimal  As, 
therefore  multiple  models. 


Table  1:  Numerical  Evaluation  of  24  x  69  Lasso-optimal  design 


9 

Important  Factors 

SVOE  Lasso  Variable  Selection 

5 

4,5,32,50,56 

(56,50,4,32,5) 

10 

3,4,6,17,23,25,26,31,37,63 

(37, 4, 31), (37,4, 31, 3), (37, 4, 31, 3, 59, 25, 17, 6) 

15 

5,8,9,20,25,26,29,33,37,44,45,46,47,60,63 

(3,37,39,46,33,5,52,24,48,58,26) 

20 

2,8,9,12,21,26,32,34,39,40, 

46,47,52,53,55,57,58,62,64,65 

(39,38,52,26,58) 

Ideally,  the  SVOE  Lasso  method  would  select  all  of  the  important  factors  and  none  of  the  unimportant 
ones.  It  does  this  for  our  q  =  5  example.  As  q  increases,  the  number  of  selected  factors  remains  roughly 
constant,  which  means  there  is  a  decrease  in  the  proportion  of  correctly  identified  factors.  Some  unimportant 
factors  also  show  up,  although  false  positives  arc  less  of  a  concern  in  screening  experiments  because  they 
can  be  eliminated  in  subsequent  stages.  Clearly,  a  much  more  extensive  empirical  investigation  would  be 
needed  to  provide  a  true  picture  of  the  design’s  screening  capabilities,  both  in  terms  of  false  positives  and 
false  negatives.  Nonetheless,  it  is  promising  to  see  that  when  the  number  of  important  factors  is  small  or 
moderate  (q  =  5  or  10),  the  SVOE  Lasso  method  selects  a  majority  or  all  of  the  important  factors  with  less 
than  half  number  of  runs  required  by  saturated  design.  This  implies  that  for  practical  screening  problems 
where  a  small  set  of  factors  can  drive  the  response,  SVOE  Lasso  may  work  well. 

With  these  promising  (though  limited)  results  giving  us  some  confidence  of  the  design,  we  ran  1000 
replications  to  obtain  sample  proportions  Y,  for  each  of  the  24  design  points.  This  took  over  5.25  days 
of  CPU  time  on  our  computing  cluster.  We  found  that  seven  of  the  24  design  points  aborted  before  all 
replications  were  complete.  After  defining  a  new  response  Y  =  I(all  reps  complete),  we  quickly  identified 
which  factor  was  causing  the  problem.  After  confirming  our  findings  by  running  MANA  in  its  GUI  mode, 
we  removed  this  factor  from  our  list  and  substituted  another.  This  points  out  a  clear  benefit  of  using  a 
screening  procedure.  We  were  able  to  identify  a  bug  in  just  5.25  CPU  days,  rather  than  the  112  CPU  days 
required  for  1000  replications  of  our  512  design  point  NOAB,  or  the  4.9  CPU  years  required  for  1000 
replications  of  an  8196  design  point  resolution  V  fractional  factorial  (Sanchez  and  Sanchez  2005). 

Our  second  screening  experiment  ran  to  completion.  There  were  two  potential  choices  for  Xopt.  The 
smallest  selects  two  of  the  first  35  factors  and  three  of  the  second  30,  the  largest  selects  six  of  the  first  35 
factors  and  five  of  the  second  30.  The  Lasso  path  for  40  <  A  <  220  appears  in  Figure  3. 
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Figure  3:  Lasso  path  for  the  24  x  65  design,  with  the  number  of  selected  terms  shown  for  each  region. 

These  results  were  somewhat  surprising,  since  we  had  anticipated  that  all  the  important  factors  would 
be  among  the  first  35.  Further  experimentation  is  underway  to  determine  an  appropriate  baseline  for  the 
study.  Whatever  these  results,  they  will  be  beneficial  for  several  stakeholders.  As  mentioned  earlier,  the 
MANA  developers  have  been  adding  features  to  expand  its  capabilities,  and  they  are  always  interested  in 
results  of  large-scale  beta  tests  that  might  indicate  either  model  bugs  (“artifacts”),  or  else  modeling  concepts 
that  would  benefit  from  additional  clarification  in  the  documentation.  Those  creating  scenarios  can  use 
screening  experiments  early  in  the  model-building  process  to  aid  in  model  verification  and  help  set  suitable 
factor  ranges.  The  international  consortium  looking  at  the  model-based  ship  design  project  will  benefit 
from  having  information  about  the  robustness  of  the  findings — we  advocate  that  it  is  better  to  experiment 
on  a  large  number  of  factors  and  seek  a  simplified  model  than  to  limit  the  factors  from  the  outset  (Sanchez 
et  al.  2012).  Finally,  by  examining  simulation  models  that  are  being  used  for  decision-making,  we  will 
continue  to  identify  challenges  that  would  benefit  from  further  methodological  research. 

4  DISCUSSION 

In  this  paper,  we  provided  a  brief  overview  of  a  newly  proposed  Lasso  method  for  factor  screening,  and  a 
new  search  criterion  for  constructing  Lasso-optimal  supersaturated  designs.  The  new  method  specifies  one 
or  more  optimal  tuning  parameters  for  Lasso,  which  determines  the  model  selection  result.  The  criterion  is 
based  on  the  probability  of  correct  variable  selection  of  Lasso.  Compared  to  the  past  work  in  supersaturated 
designs,  our  work  uses  the  analysis  method  to  direct  the  experimental  design  stage,  and  the  preliminary 
numerical  evaluation  shows  very  promising  results.  The  method  is  specifically  intended  for  the  initial  study 
of  complex  systems  with  many  potential  factors. 

The  example  in  the  paper  involves  an  agent-based  simulation  of  maritime  escort  operation.  The  insights 
we  gain  in  this  study  may  continue  to  aid  in  model  verification  and  validation  efforts — we  have  already 
identified  one  model  bug  by  using  the  small  Lasso-optimal  screening  design.  The  broader  implications  are 
also  important:  the  operational  insights  from  Kaymal’s  investigation  (and  related  studies)  may,  in  the  near 
future,  determine  how  nations  decide  to  structure  their  naval  forces. 

In  summary,  preliminary  results  on  the  supersaturated  designs  are  promising,  both  for  artificially 
generated  response  surfaces,  and  for  a  more  complex,  agent-based  simulation.  Supersaturated  screening 
designs  may  be  of  even  greater  interest  for  larger  models:  many  simulations  addressing  defense  or  homeland 
security  issues  involve  thousands  of  potential  factors. 
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